%% 清理
clear; clc; close all;

%% 1. 定义初始条件函数 eta(x)
eta = @(x) exp(-20*(x-2).^2) + exp(-(x-5).^2);

%% 2. 精确解函数
exact_solution = @(x, T) eta(x - T);

%% 4. 设置参数
h = 0.05;
T = 17;
x_min = -30; x_max = 50;
x_big = (x_min:h:x_max)';
mu_list = [0.8, 1.0];

schemes = { @leapfrog, @upwind, @beam_warming, @lax_friedrichs, @lax_wendroff };

plot_cfg = {
    {@leapfrog, 0.8, 'Leapfrog (μ=0.8)'}, ...
    {@upwind, 0.8, 'Upwind (μ=0.8)'}, ...
    {@beam_warming, 0.8, 'Beam‑Warming (μ=0.8)'}, ...
    {@lax_friedrichs, 0.8, 'Lax‑Friedrichs (μ=0.8)'}, ...
    {@lax_wendroff, 0.8, 'Lax‑Wendroff (μ=0.8)'}, ...
    {@lax_wendroff, 1.0, 'Lax‑Wendroff (μ=1.0)'}, ...
    {@leapfrog, 1.0, 'Leapfrog (μ=1.0)'}};


figure('Position',[100 50 800 1200]);
subplot(4,2,1);
x_init = -2:h:8;
plot(x_init, eta(x_init),'k-','LineWidth',1.5);
title('Initial Condition'); legend('Initial');
grid on;

cnt = 2;
for cfg = plot_cfg
    method = cfg{1}{1};
    mu = cfg{1}{2};
    label = cfg{1}{3};
    
    k = mu * h;
    nsteps = floor(T / k);
    u0 = eta(x_big);
    u_num = method(u0, mu, nsteps);
    u_ex = exact_solution(x_big, T);
    
    idx = x_big>=15 & x_big<=25;
    x_plot    = x_big(idx);
    u_num_p   = u_num(idx);
    u_ex_p    = u_ex(idx);
    
    subplot(4,2,cnt);
    plot(x_plot, u_ex_p, 'k-','LineWidth',1.5); hold on;
    plot(x_plot, u_num_p,'r-','LineWidth',0.8);
    title(label);
    
    grid on;
    cnt = cnt+1;
end

suptitle('Reproduce Figure 12.8: Accuracy of Various Methods for the Advection Equation');


function u_new = leapfrog(u0, mu, nsteps)
    N = numel(u0);
    u_old = u0;
    u_new = u0;
    u_new(2:N-1) = u0(2:N-1) - mu*(u0(2:N-1) - u0(1:N-2));
    for n = 2:nsteps
        u_temp = u_new;
        u_new(2:N-1) = u_old(2:N-1) - mu*(u_temp(3:N) - u_temp(1:N-2));
        u_old = u_temp;
    end
end

function u_new = upwind(u0, mu, nsteps)
    u_old = u0;
    N = numel(u0);
    u_new = u_old;
    for n = 1:nsteps
        u_new(2:N) = u_old(2:N) - mu*(u_old(2:N) - u_old(1:N-1));
        u_new(1) = u_old(1);
        u_old = u_new;
    end
end

function u_new = beam_warming(u0, mu, nsteps)
    u_old = u0;
    N = numel(u0);
    u_new = u_old;
    for n = 1:nsteps
        u_new(3:N) = u_old(3:N) ...
            - mu/2*(3*u_old(3:N) - 4*u_old(2:N-1) + u_old(1:N-2)) ...
            + mu^2/2*(u_old(3:N) - 2*u_old(2:N-1) + u_old(1:N-2));
        u_new(1:2) = u_old(1:2);
        u_old = u_new;
    end
end

function u_new = lax_friedrichs(u0, mu, nsteps)
    u_old = u0;
    N = numel(u0);
    u_new = u_old;
    for n = 1:nsteps
        u_new(2:N-1) = 0.5*(u_old(3:N) + u_old(1:N-2)) - mu/2*(u_old(3:N) - u_old(1:N-2));
        u_new([1, N]) = u_old([1, N]);
        u_old = u_new;
    end
end

function u_new = lax_wendroff(u0, mu, nsteps)
    u_old = u0;
    N = numel(u0);
    u_new = u_old;
    for n = 1:nsteps
        u_new(2:N-1) = u_old(2:N-1) ...
            - mu/2*(u_old(3:N) - u_old(1:N-2)) ...
            + mu^2/2*(u_old(3:N) - 2*u_old(2:N-1) + u_old(1:N-2));
        u_new([1, N]) = u_old([1, N]);
        u_old = u_new;
    end
end

